import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 定义目标函数（均方误差）
def objective_function(p, actual_defective_rate):
    return np.mean((p - actual_defective_rate) ** 2)

# 粒子群优化算法
class PSO:
    def __init__(self, n_particles, n_dimensions, objective_func, actual_defective_rate, bounds, max_iter=100, w=0.5, c1=1, c2=2):
        self.n_particles = n_particles
        self.n_dimensions = n_dimensions
        self.objective_func = objective_func
        self.actual_defective_rate = actual_defective_rate
        self.bounds = bounds
        self.max_iter = max_iter
        self.w = w
        self.c1 = c1
        self.c2 = c2
        self.particles = np.random.uniform(bounds[0], bounds[1], (n_particles, n_dimensions))
        self.velocities = np.zeros((n_particles, n_dimensions))
        self.pbest_positions = self.particles.copy()
        self.pbest_values = np.array([objective_func(p, actual_defective_rate) for p in self.particles])
        self.gbest_index = np.argmin(self.pbest_values)
        self.gbest_position = self.pbest_positions[self.gbest_index]
        self.gbest_value = self.pbest_values[self.gbest_index]

    def update_particles(self):
        for i in range(self.n_particles):
            r1, r2 = np.random.rand(2)
            self.velocities[i] = (self.w * self.velocities[i] +
                                  self.c1 * r1 * (self.pbest_positions[i] - self.particles[i]) +
                                  self.c2 * r2 * (self.gbest_position - self.particles[i]))
            self.particles[i] += self.velocities[i]
            self.particles[i] = np.clip(self.particles[i], self.bounds[0], self.bounds[1])
            fitness = self.objective_func(self.particles[i], self.actual_defective_rate)
            if fitness < self.pbest_values[i]:
                self.pbest_values[i] = fitness
                self.pbest_positions[i] = self.particles[i]
            if fitness < self.gbest_value:
                self.gbest_value = fitness
                self.gbest_position = self.particles[i]

    def run(self):
        gbest_values = []
        for _ in range(self.max_iter):
            self.update_particles()
            gbest_values.append(self.gbest_value)
        return self.gbest_position, self.gbest_value, gbest_values

for i in range(12):
# 实际次品率
    actual_defective_rate = 0.1#0.05， 0.1， 0.2

    # 粒子群参数
    n_particles = 1000
    n_dimensions = 1
    bounds = [0, 1]
    max_iter = 100

    # 初始化并运行PSO
    pso = PSO(n_particles, n_dimensions, objective_function, actual_defective_rate, bounds, max_iter)
    best_position, best_value, gbest_values = pso.run()

    print(f"估计的次品率为：{best_position[0]}")
    print(f"目标函数值（MSE）为：{best_value}")

    # 抽样检测方案
    sample_size = 98  # 根据之前的计算，样本大小为98

    # 模拟抽样检测
    samples = np.random.choice([0, 1], size=sample_size, p=[1 - best_position[0], best_position[0]])
    defective_count = np.sum(samples)
    defective_rate = defective_count / sample_size

    print(f"抽样检测结果：次品率为 {defective_rate}")

    # 决策规则
    if defective_rate > 0.15: # 0.1，0.15，0.25
        print("在95%的信度下，次品率超过标称值，拒收这批零配件。")
    else:
        print("在95%的信度下，次品率未超过标称值，接收这批零配件。")

# # 可视化
# plt.figure(figsize=(10, 6))
# plt.plot(gbest_values, label='Global Best Value')
# plt.xlabel('Iteration')
# plt.ylabel('MSE')
# plt.title('PSO Optimization Process')
# plt.legend()
# plt.grid(True)
# plt.show()
